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ABSTRACT 

o : 

i I present a model-independent spherically symmetric density estimator to be used in the 

cross-correlation of imaging catalogs with objects of known redshift. The estimator is a simple 
O . modification of the usual projected density estimator, with weightings that produce a spherical 

| aperture rather than a cylindrical one. 

C*") , Subject headings: large-scale structure of the universe — methods: statistical 

> : 

1. Introduction 

00 . 

. Measuring galaxy properties as a function of their local environment is a central task of observational 

t-H ' extragalactic astronomy. Doing so requires a measure of density. With redshift surveys, one can estimate 

densities by counting other spectroscopic galaxies in a redshift-space window around each primary object. 
However, imaging catalogs are generally much deeper than spectroscopic catalogs, which suggests that cross- 
correlating the imaging catalog with the spectroscopic catalog should yield useful densities even for the 
faintest spectroscopic objects. 

o : 

The presence of at least one spectroscopic redshift in each pair is of great utility because it means 
that one can immediately map angular separations to physical distances and derive intrinsic properties 
(i.e. luminosities) for both the spectroscopic object and the correlated objects from the imaging catalog. 
, ^ | This opportunity has been used extensively in the early study of galaxy clustering (Davis et al. 1978; Yee & 

^ ■ Green 1987; Lilje & Efstathiou 1988; Saunders et al. 1992) and in the study of dwarf galaxies in groups (e.g., 

Ferguson & Sandage 1991) and in the field (Phillipps & Shanks 1987; Vader & Sandage 1991; Lorrimer et 
al. 1994; Loveday 1997). Of course, the correlated objects are still seen in projection, with a mix of spatial 
separations at each transverse separation. 

Angular correlations can be inverted to spatial correlations, and hence to a form of density, by the 
assumption of isotropy, and many of the above studies did this by assuming that the correlation functions 
are power laws in scale (e.g., Phillipps 1985). Saunders et al. (1992) and Loveday (1997) go one step further 
by using the assumed power-law as an optimal filter. Fall & Tremaine (1977) suggest that the general 
inversion can be done with a smoothing filter. Baugh & Efstathiou (1993) use a regularized Lucy's iteration 
to do a similar inversion, while Dodelson & Gaztahaga (2000) and Eisenstein & Zaldarriaga (2001) use other 
smoothing priors. 

In the study of galaxy environmental dependences on small scales, it is useful to pursue a more model- 
independent measure of the density. On scales of 1 Mpc, it is quite possible that the correlation functions 
will deviate from power laws (and certainly from a uniform power law) in ways that depend on luminosity, 
star- format ion rate, or other variables (e.g., Scranton 2002). 
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Hcre I describe a method to recover density estimates in spherically symmetric real-space windows from 
the cross-correlation of imaging and spectroscopic catalogs independently of the shape of the correlation 
function. The resulting formulae are a simple alteration of the usual background subtraction methods. The 
method can be subdivided so as to yield noisy density estimates for individual spectroscopic objects that 
can then be averaged according to whatever subclasses one might desire. 

2. Angular Cross-correlations 
2.1. Definitions 

We begin with two samples of objects, one with spectroscopic redshifts and the other without. We bin 
the spectroscopic sample into thin redshift bins (indeed, we can consider each spectroscopic object separately) 
and select a subsample of the imaging catalog using the known redshift. For example, one might select a 
sample in a particular luminosity range (i.e. the objects would have the desired luminosity were they at 
the spectroscopic redshift). We adopt a flat-sky coordinate system in which the transverse directions arc 
measured as distance R at the given redshift and the linear direction is measured as distance Z where Z = 
is at the redshift. We write the three-dimensional position (R, Z) as r. 

Of course, the imaging subsample samples a range of Z, not just Z — 0. We describe the homogeneous 
density of the sample as (fi(Z), where this is the number per unit Z and per unit transverse area. The areal 
number density is 

/oc 
dz 4>{Z). (l) 
-oo 

If the true spatial cross-correlation between the spectroscopic object and the imaging subsample is £i s (r), 
then the angular correlation is 

Wis(R) = \ r dZ 0(Z)UVZ? + R 2 ) « ^ j X dZ UVZ 2 + R 2 ) (2) 

where the last equality assumes that the selection function <fr is essentially constant over the scale on which 
£i s is not negligible. We also assume that the angular diameter distance is constant in this region. We denote 
</>(0) as 4>o- Equation (2) is a simple form of Limber's equation (Limber 1953; Groth & Peebles 1977). 



2.2. Deprojection 

It is well-known that equation (2) can be inverted as an Abel integral (von Zeipcl 1908), so that 

, n dR dw is (R) (3) 

However, because this involves a derivative of the measured w; s (i?), it is noisy. 
We address this by measuring only an integral of £i s 

1 r°° 

A = - Anr 2 dr £ is {r)W(r) (4) 
* Jo 

where V — Ai:r 2 dr W{r). W(r) is our smoothing window. A has a very useful physical meaning: it 
is the average overdensity of objects from the imaging catalog in the neighborhood (as defined by W) of a 
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spectroscopic object. Note that A is defined in 3-space with a spherically symmetric window; it is not a 
projected quantity. 

Inserting (3) into (4) and switching the limits of integration yields 

0o Vir J dR J Vi? 2 - r 2 <j> V J dR y ' K > 

where 

For bounded W(r), F(0) = and F(R) — » F/27r 2 i? for large i?. Integrating by parts yields 

A = - — / 2irRdRw is (R)G(R) (7) 
<PoV J 

G <*> s is- < 8 > 

A constant w is integrates to A = 0, so adding a constant to W[ s doesn't change A. G(R) — ► — y/27r 2 i? 3 for 
large i?. 

Now let us consider our measurement of w- ls (R). In a small radial bin from R to R + dR, we would 
estimate 1 + Wi s (R) as the ratio of the observed counts of pairs in that radial range to the expected number. 
The expected number is 2irnRdR. Treating the integral as a Riemann sum in which the bins are so small 
as to contain or 1 observed pair leads to the conclusion that 

* = £ £ ^ E o» 

where the sums are over the objects in spectroscopic subsample and the imaging subsample, respectively, 
N sp is the number of spectroscopic objects, and Rjk is the transverse separation of the j th spectroscopic 
object to the k th imaging object. 

An important point is that we can now treat each spectroscopic object separately, yielding the following 
noisy measure of the overdensity around object j: 

k£{im} 

Note how simple this formula is: one counts the imaging objects, weighting by G(R), and divides by the 
expected number of objects in the real-space window (V<po). We can recover the average density around any 
subset of the spectroscopic sample simply by averaging the selected Aj . 

It is interesting to compare equation (10) to a more conventional background-subtraction method in 
which one sums all of the objects in an angular aperture and subtracts an appropriately scaled value of the 
areal density averaged over the entire survey. This would correspond to a G function that was constant 
and positive for R less than the aperture and then constant and negative for all greater R. The difference 
is that this background subtraction would give an estimate for the density in a cylindrical region, in which 
the axis of the cylinder lies along the line of sight and is much longer than the radius of the cylinder. The 
resulting density would sample the correlation function £; s at a wide range of radii. The formula given here 
creates a compact region in all three dimensions. Some workers (e.g., Gaidos 1997; Valotto et al. 1997) have 
used annular regions for the determination of the background. This truncates the cylindrical region in some 
fashion, but the detailed effects were not assessed. 
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2.3. Gaussian Windows 

For a useful and illustrative example, we will treat the case in which the window is assumed to be a 
Gaussian Wir) = cxp(-r 2 /2a 2 ). The volume of the window is V — (27r) 3 / 2 a 3 . Then 

2 f R r 2 e -r 2 /2a 2 n2 

nR) = -J o d r 7 == f = T e-s [ I ( S )-I 1 ( S ) ] (11) 

where s = i? 2 /4a 2 and I n are modified Bessel functions of the first kind (Gradshteyn & Ryzhik 2000, 3.364.1). 
Then 

G(R) = e- s [I (s) - 2sl (s) + 2sh{s)] . (12) 

The weighting function G is shown in Figure 1. The function is smooth, with a positive peak at R s=s a 
and a broader negative peak at R w 3a. As expected, to estimate the density averaged over the window W, 
one is counting the nearby objects and subtracting a count of objects slightly further away. 



2.4. Boundaries and Masks 

Equation (7) requires an integration over all radii and since G cx R~ 3 at large radius, this integration 
converges as ~ 1/R. If one splits the integral in equation (7) at some radius i? max , inside of which one will 
do the counting in equation (9), then the residual error from larger radii is 



/>OC 

/ 2nRdR [1 + w is (R)]G(R) (13) 



Here, we have to include the constant background term because it will no longer integrate to zero. We 
can estimate the correlated term by treating Wi s as a power-law w ls (R) = Wi s (R maK )(R max / R) a , where 
Wi s (i? max ) is the value at the truncation radius and a is the slope of the power-law, typically 0.7-0.8 in past 
measurements. Then, using G(R) w — V/2n 2 R 3 in the second term, we find 



2^F(i? max ) w is {R max ) 1 



(14) 



V 7ri? max 1 + q_ 

Taking £; s cx r _1_a , we have 

^ Wis{ R) = R^R) 1 ^! (15) 



One finds that A res = — 2nnF(R mayi ) / 4>qV — 0.70£i s (i? max ) for a = 0.75. Applying this as a correction does 
introduce some model dependence on the form of w- m , but one can pick i? max so that the correction is rather 
small. In this way, one can avoid summing over all pairs of spectroscopic and imaging objects. 

The circular region R < i? max may still include regions that are outside the survey or its mask. If 
one writes $(i?) as the fraction of the annulus of radius R that is within the survey, then one can weight 
the counts in equation (9) by 1/$(_R). However, $(i?) may be expensive to compute for each spectroscopic 
galaxy. An alternative method is to generate a catalog of random points that are uniformly distributed 
outside of the survey region and then add to equation (10) the sum over those points closer than i? max , 
weighting by G(R)[1 + w ls (R)} and the ratio of n to the random catalog surface density. This effectively 
interpolates over the masked regions. If one treats W\ S as a function of fixed shape but unknown amplitude, 
then one can do the sum of the clustered term separately and renormalize after one has determined the 
amplitude from the co-added Aj. Of course, since one is in some fashion assuming the answer, one should 
exclude spectroscopic galaxies for which the masked contribution is large. 
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2.5. Errors and Optimal Weighting 

When averaging the individual estimates Aj from a sample of spectroscopic objects, one can achieve 
more precision by using non-uniform weighting in the mean. The optimal weighting is inverse variance; 
however, one should remember that approximations to inverse variance may be much simpler to compute 
and yet still close to optimal. One should never weight by variances that are derived from single data points, 
e.g. the square root of the observed number of companions to a given galaxy, because this will bias the 
resulting mean. 

The statistical variance of the Aj estimator comes from two sources, the clustering of the galaxies and 
shot noise. The clustering term involves the three-point correlation function, with one galaxy from the 
spectroscopic sample and two from the imaging sample. If the first of these is at the origin and the latter 
two are at r\ and r 2 , then we denote the three-point correlation function as ((ri, r 2 ). The expected densities 
at two points given a spectroscopic galaxy at the origin is (Peebles 1980) 

(p(?i)p(f2)> - <t>(Zi)<f>(Z 2 ) [1 + Un) + Ur 2 ) + - f 2 |) + C(ri,r 2 )] (16) 

where £; s is the two-point correlation between imaging and spectroscopic galaxies and £u is the correlation 
of two imaging galaxies. 



The clustering contribution to the variance of Aj is then 

Ki(Aj) = J d 2 Ri J d 2 R 2 G(R\)G(R 2 ) J dZidZ 2 



x 



<HW(z 2 ) 



- r 2 |) + &(ri)&(r 2 )] ( 17 ) 



V *V* 



The last term is simply (Aj) 2 . The terms in £; s in equation (16) integrate to zero. The £ term involves 
all three objects, so at the level of approximation in equation (2), we can assume <fi(Zi) = 4>{Z 2 ) = 4> - 
However, the term represents correlated imaging galaxies that are uncorrelated with the spectroscopic 
object, so they may have Z far from and <j>{Z) ^ cfio- Integrating over Z will yield the angular correlation 
of the imaging catalog. Indeed, for Z ^ 0, the angular separations implicit in our definition of the transverse 
coordinate R correspond to different physical scales. Doing this correctly again yields the angular correlation 
function. We thus simplify equation (17) to 

V cl {Aj) = V 2pt + V 3pt (18) 

V 2pt = [^^j j d 2 R 1 j d 2 R 2 G(R 1 )G(R 2 )w u (\R 1 -R 2 \) 

Vz P t = J^I^G(R 1 )G(R 2 )[an,f 2 )-Uri)Ur2)} 

The first term can be done quickly by Fourier methods, as can the second if one adopts the hierarchical 
ansatz (Groth & Peebles 1977) to write C(ri,r 2 ) = Q{6s(ri)£ is (r 2 ) + [£ is (n) + £ is (r- 2 )]&i(|fi - r 2 )|)}. The 
two-dimensional Fourier transform of G(R) is 

J d 2 Re lin G(R) =4 drrsm(kr)W(r) (19) 

Note that while V 3pt is a contribution to the variance about the mean Aj, it is not necessarily noise! Much 
of it is the density around the particular object, which of course has scatter from the mean. 



- 6 - 



The shot noise or Poisson contribution to the variance is based on the expected counts, including 
clustering, which are n[l + Wi S (R)]. The variance is 

2 

n[l + w is {R)]. (20) 

We refer to these two terms as the homogeneous and clustered shot noise, respectively. 

For the Gaussian window (eq. [12]), the homogenous shot noise becomes 2a 2 n/ (V<fio) 2 . Figure 1 shows 
the contribution per radial bin to the variance in the shot noise. Essentially all of the shot noise arises at 
R < 1.5a; in other words, the fact that one is subtracting the background with a region at moderate radius 
rather than the entire sample adds little extra noise to the density estimator. 

The four above contributions — V2pt, V3 pt , homogeneous shot noise, and clustered shot noise — all have 
different scalings with the depth of the survey, the size of the window, and the clustering strength. For 
a typical survey thickness L = n/(f>o and a typical window radius a, the contributions to the variance are 
roughly (i/20a)A, A 2 , (L/125a)(l/a 3 </>o), and (A/32)(l/a 3 ^>o), respectively. The numerical coefficients are 
for illustration only 1 . The 2-point clustering term dominates on large scales; the clustered shot-noise on 
small scales. 

The process of averaging many Aj into a mean A introduces additional error terms from the correlations 
of the spectroscopic galaxy positions, including contributions from the four-point correlation function. This is 
not surprising because A is an integral of Wi s (R), whose covariance normally involves the four-point function. 
Neglecting these new terms in favor of the equations above corresponds to the assumption that the shot noise 
of the spectroscopic sample exceeds its clustering (see the expansions of Bernstein 1994; Hamilton 1997). 
This is often a good assumption on small scales, particularly if one is considering only a small subset of the 
spectroscopic sample. 

A crucial assumption of the analysis is that the objects that are uncorrelated but in close projection 
with the spectroscopic object are statistically identical to those in other parts of the sky. Magnification from 
weak lensing can violate this assumption in principle (e.g., Valotto et al. 1997). Moreover, selection biases 
in cluster catalogs owing to superposition of unrelated structures (e.g., Valotto et al. 2001) are not reduced 
by this method. 



y sn (A,) = / d 2 R 



G(R) 



2.6. From Kernels to Windows 

If we are given G(R), we can use the Abel integral to find the corresponding W(r): 

W(r) = - f d R ^B= . (21) 
r J Vr 2 - R 2 

If W(r) is to be bounded and V ^ 0, G(R) must have an asymptotic form of — V/2n 2 R 3 . In particular, 
this means that one cannot have G = for all large radii, thereby avoiding summations over large pair 
separations. 



lr Fhe numerical factors are computed for the Gaussian window with the assumptions of a power-law correlation function 
scaling as r -18 with uniform bias for all galaxies, an a = —0.9 Schechter luminosity function with imaging catalog luminosity 
cuts between 0.4L* and 2.5L* , the hierarchical form of f with Q = 1.3, and a non-expanding cosmology (Euclidean metric, no 
K-corrections). 
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2.7. Other Windows 



As one might expect, the spherical tophat is easily calculated but poorly behaved. Taking W(r) = 1 
for r < b, we find 

1 R < b 



G(R) 



■ -i b 
sin — — 



R> b 



(22) 



R VR 2 -b 2 _ 

This diverges as R — > 6 from above. The singularity is integrable in A but divergent in the shot noise, 
making the tophat a poor choice. 

Another simple choice that retains the advantage that W = for r > b is W(r) = 1 — r 2 /b 2 for r < b. 
We find 



F(R)=l { 



R 2 


3R 4 


2 




2 




7T 





r 2 m A 

2 8b 2 



sin 



and 



G(i2) = < 



3i? 2 
26 2 



3Vfl 2 - b 2 
26 



3_R^ 

26 2 



i? 



R 

R<b 
R>b 



R<b 



R>b 



(23) 



(24) 



This is continuous at R = b but has somewhat worse shot noise properties than the Gaussian window. 



3. Conclusions 

I have presented a model-independent estimator for the spherically averaged overdensity of imaging 
catalog objects around spectroscopic objects. The method is simple to apply; one can view it as an alteration 
of standard background subtraction methods to yield spherical apertures rather than cylindrical ones. The 
method does not require an explicit inversion of the angular cross-correlation function to a spatial correlation 
function, although clearly if one finds oneself measuring the densities in multiple apertures of different sizes, 
one is effectively reverting back to an inversion method. The primary advantages of the new method compared 
to integrating the output of an inversion (i.e., computing eq. [3] and then eq. [4]) are that one does not need 
a smoothing prior and that the statistic can be applied to each spectroscopic galaxy independently (eq. [10]), 
leaving one free to sum the results post facto across as many spectroscopic subsamples as one desires. Like 
other angular methods, the new density estimator is unaffected by redshift distortions, which gives it an 
advantage on small scales over density estimation from spectroscopic catalogs. 

Past work has assumed power-law correlation functions to recover spatial densities. This allows one 
to achieve higher signal-to-noise ratio and hence is a better choice for some applications. However, when 
probing the dependence of galaxy properties on small-scale environment, the model independence of the 
density estimator presented here is a valuable advantage. With today's large surveys (e.g., York et al. 2000; 
Colless et al. 2001), statistical precision is sometimes less precious than systematic control. 

While it is clear that one can profitably consider the dependence of average density on the properties of 
the spectroscopic objects (Hogg et al. 2002), it is worth pointing out that one can also derive densities for 
subdivisions of the imaging catalog. For example, one can probe the color and/or luminosity distribution of 
objects within a spherical aperture of a particular set of spectroscopic objects (e.g., galaxies or clusters). A 
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speculative application would be to couple this approach to galaxy-galaxy weak lensing mass estimates. If 
one had a mass estimate for each imaging object, assuming the spectroscopic redshift, then one could find 
the masses of objects that are correlated with the spectroscopic tracer. 

I thank David Hogg, Jon Loveday, Tim McKay, Ann Zabludoff, and Dennis Zaritsky for useful discus- 
sions. D.J.E. was supported by National Science Foundation (NSF) grant AST-0098577 and by a Alfred P. 
Sloan Research Fellowship. 
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Fig. 1. — (solid line) The integration kernel G(R) for the Gaussian window, (dashed line) RG(R), which 
is the weight per radial bin for a uniform background, (dot-dashed line) RG 2 (R), which is the shot noise 
variance per radial bin for a uniform distribution. 



